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The dynamics of dissipative soft-sphere gases obeys Newton's equation of motion which are com- 
monly solved numerically by (force-based) Molecular Dynamics schemes. With the assumption of 
instantaneous, pairwise collisions, the simulation can be accelerated considerably using event-driven 
Molecular Dynamics, where the coefficient of restitution is derived from the interaction force be- 
tween particles. Recently it was shown, however, that this approach may fail dramatically, that is, 
the obtained trajectories deviate significantly from the ones predicted by Newton's equations. In 
this paper, we generalize the concept of the coefficient of restitution and derive a numerical scheme 
which, in the case of dilute systems and frictionless interaction, allows us to perform highly effi- 
cient event-driven Molecular Dynamics simulations even for non-instantaneous collisions. We show 
that the particle trajectories predicted by the new scheme agree perfectly with the corresponding 
(force-based) Molecular Dynamics, except for a short transient period whose duration corresponds 
to the duration of the contact. Thus, the new algorithm solves Newton's equations of motion like 
force-based MD while preserving the advantages of event-driven simulations. 

PACS numbers: 45.50.Tn, 45.70.-n 



I. INTRODUCTION 

Modelling granular systems of frictionless spheres 
branches into two fundamental different approaches: 
Hard- and soft-sphere models. The dynamics of soft 
spheres are governed by the pairwise interaction forces 
between contacting particles as a function of the rela- 
tive particle positions and velocities as well as material 

parameters, Fij = Fij (^fi^fj^fi^fj^. The dynamics of 

a many-particle system is then obtained by numerically 
solving Newton's equation of motion for all degrees of 
freedom, which was termed Molecular Dynamics (MD), 
e.g. [T. The first MD simulations of granular systems 
(in the engineering literature also called Discrete Element 
Method - DEM) range back to pioneering work by Cun- 
dall, Walton, Haff and others, e.g. i2H5j. An overview 
of the force models specific for granular particles can be 
found in [MS]. 

In contrast to soft-sphere models, in hard- sphere mod- 
els the collisions are assumed to occur instantaneously 
which allows to consider the dynamics of hard sphere 
systems as a sequence of independent binary collisions. 
Except for collisions where the velocities change instan- 
taneously, the particles follow ballistic trajectories, possi- 
bly under the influence of external fields like gravity. The 
hard sphere model is the foundation of both, Kinetic The- 
ory of granular matter based on the Boltzmann equation 
e.g. ^9 -11 , and event-driven Molecular Dynamics (eMD) 
of granular matter, e.g. [ 121414] . 

The collision of two hard spheres of velocities fi and 
fj implies an instantaneous exchange of momentum: 

(ri' - f;') • e; = -s„ (f;° - r^) ■ e ° (i) 

with the time dependent inter center unit vector eV = 
{vi — Vj) I \Ti — rj\ and the coefficient of normal restitu- 



tion En- Upper index denotes values just before the 
collision, primed values denote post-collisional values. 
Unlike the velocities, the particles' positions remain un- 
changed because of the instantaneous character of the 
collision, therefore, 

K ^ e ° (2) 

and Eq. ([T]) reduces to 

(?i - f;') . e ° = -£„ - -rf) ■ . (3) 

Equation ([3| relating the pre- and post-collisional veloc- 
ities is the governing equation of eMD. Given a certain 
granular system may be described by the hard-sphere 
model, eMD allows for a vast increase of numerical effi- 
ciency as compared with corresponding MD simulations. 
For a very efficient implementation of eMD see [15^ . 

Despite of eMD's great numerical performance, the 
hard-sphere model is a simplification of physical reality: 
instantaneous changes of velocity imply infinite delta- 
shaped forces while forces between colliding physical ob- 
jects are always finite which implies finite contact du- 
ration. Therefore, the applicability of the hard-sphere 
model for eMD simulations of granular systems must 
be checked. One obvious precondition for eMD is low 
enough particle number density such that the frequency 
of three-particle contacts can be neglected as compared 
to the frequency of pair collisions. Obviously, this is not 
given for slow flows with long-lasting contacts. 

A natural way to check the validity of the hard-sphere 
approximation in the dilute limit is the following: The 
coefficient of normal restitution as a function of mate- 
rial parameters and relative impact velocity may be ob- 
tained from analytically integrating Newton's equation 
of motion for the central collision of an isolated pair of 
particles using the known interaction force which is also 
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a function of material properties and impact velocity, e.g. 
p^BHTO] . Performing now MD simulations using the inter- 
action force and eMD simulations using the correspond- 
ing expression for the coefficient of restitution, one may 
expect identical trajectories. However, recently it was 
found that these trajectories may deviate significantly 
for a vast range of materials, collision geometries and 
impact velocities [20, 21 , in particular for oblique im- 
pacts which concerns the majority of impact geometries 
[20] in a Molecular Chaos situation. Consequently, even 
for dilute systems, the hard-sphere approximation may 
fail dramatically. This effect may be attributed to the fi- 
nite duration of collisions in physical systems which does 
not allow for the assumption Eq. ([2|. 

Consequently, on one hand we have the stunning effi- 
ciency of eMD based on the hard-sphere model. On the 
other hand there is the universality and physical correct- 
ness of the soft-sphere model leading to MD. Combining 
the advantages of both approaches is a highly desired 
aim. Concerning simulation techniques an attempt is to 
discretize the (smooth) interaction potentials. As this 
idea was originally developed for liquids [22 , 23 recently 
it was also applied to granular systems [24-26 . On the 
theoretical side there are perturbation theories, extend- 
ing hard sphere models [27H29] . 

In this work we derive an algorithm for the event- 
driven simulation of smooth spheres which does not rely 
on Eq. ([2|. By extending the concept of the coefficient 
of restitution, we map the correct Newtonian dynamics 
of soft spheres to instantaneous events. We show that 
for dilute systems of frictionless particles the presented 
method allows for a correct computation of the trajec- 
tories (as MD) while preserving the efficiency of event- 
driven simulations. 

This simulation method applies to a wide range of 
particle interaction forces. Here we demonstrate it for 
the case of two important examples: The linear dashpot 
model and viscoelastic spheres. Unlike the original eMD 
method, we show that in these cases the trajectories ob- 
tained by eMD agree perfectly with the MD results. 



II. COLLISION OF SPHERES 

Consider two colliding spheres of masses and rrij lo- 
cated at fi{t) and fj{t) and travehng with velocities fi{t) 
and rj{t). With the interaction force F, their motion is 
described by 



m^^r = F, MR = 



where 



miTrij 
rUi + rrij 



(4) 



(5) 



are the center of mass coordinate, the relative coordinate 
and the effective mass, respectively. The center of mass 



moves due to external forces such as gravity and sepa- 
rates from the relative motion which in turn contains the 
entire collision dynamics. 

For frictionless particles, the interaction force acts in 
the direction of the inter-center unit vector, F = FnCr. 
During the collision the (orbital) angular momentum is 
conserved which allows for the definition of the constant 
unit vector cl: 

L = meff r X r = Lcl - (6) 
Thus, with the coordinate system S spanned by 

(7) 



and with its origin in the center of mass the collision 
takes place in the et^-e^-plane [46] . In the collision plane 
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FIG. 1. Illustration of the used polar coordinates (see text) 

we formulate the equation of motion in polar coordinates 
{r,(p} (see Fig.[l}: 

meff r'^(p = L, meff r = Fc ^ F^ = meff rcp'^ + F^ , (8) 

with the centrifugal force Fc. Together with the inital 
conditions 

r(0)=r^ r(0)=r^ ^(0) = , (9) 

Eq. (|8| fully describes the collision dynamics for an arbi- 
trary normal force F^. The collision terminates at time 
t = T where [H Hi] 



f(r) > and Fn = 0. 



(10) 



Measuring time in units of T, length in units of X 
and angles in units of and using the dimensionless 
quantities 



r ~ t ^ Lp 
r = — , t = 7^ and (p ■ 



(11) 



we obtain the scaled form of the equation of motion, 
Eq. (p|: 



difi 












dt2 


V dt 



(12) 
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where X and T are length and time scales typical for 
the given normal force F^, ^ is an arbitrary scale for 
measuring angles and c^^ reads 



(13) 



The corresponding dimensionless initial conditions read 
^(0)=0, r(0) = ^ and ^(0)=f(0)|. (14) 



According to Eq. ( 10 ) the scaled contact duration reads 
f = r/T. 



III. 



COEFFICIENT OF RESTITUTION VERSUS 
MATRIX OF RESTITUTION 



Solving the scaled equation of motion Eq. ( 12 ) for the 
initial conditions Eq. (14) in the time interval < t < f, 
that is, from the beginning of the collision at time t = 
until its end at t = f (see Eq. (10)) we obtain the post- 



collisional values (^(f), (^(f), f(f) and r(f), which deter- 
mine the state of the system at the end of the collision. 

Note that for the special case of central collisions with 
vanishing angular momentum, the state would be fully 
described by f(f) and f(f) as the other values vanish or 
are invariant. Together with the hard-spere assumption, 
Eq. (§, we are left with only f(f) which allow us to car- 
acterize the collision by a single number, the coefficient 
of restitution, 



f(f) 

m 



(15) 



Therefore, Sn ^ [0, 1] for central collisions. 

Equation ([l5| provides the link between the hard- 
sphere and the soft-sphere models and, correspondingly, 
between MD and eMD since it relates the coefficient of 
restitution with the specific interaction force. The ana- 
lytical solution of Eq. ( 15 ) is frequently non-trivial, even 



for rather simple forces as the viscoelastic Hertz force 
[171 [H]- As a result from the solution of Eq. (15) we 



obtain the coefficient of normal restitution as a function 
of the force's material specifics, particle sizes and impact 
rate. 

Obviously, only for the special case of central collisions 
of vanishing duration, is sufficient to characterize col- 
lisions since otherwise (^(f), (p{f) and f(f) do not vanish. 
It was shown for ordinary material and impact param- 
eters, that the mentioned post-collisional quantities are 
not negligible [20]. If one overrides this fact, the coeffi- 
cient of restitution must depend on the impact parameter 
d (see Fig. |2] below). Depending on d it can adopt even 
negative values [21]. Therefore, we believe that know- 
ing the coefficient of restitution, e^, is not sufficient to 
perform particle simulations. 

Following the previous arguments, besides the ordi- 
nary coefficient of normal restitution, we define further 



coefficients which together characterize the collision com- 
pletely. These are 



f(f) 
f(0) ' 



(16) 



which stands for distance of the colliders at the end of 
the collision. Naively one could believe = 1 since the 
particles lose contact when \ri — rj \ = Ri-\-Rj. However, 
as shown in [16[|17], the latter condition is not correct and 
leads to erroneous attractive forces even if the interaction 
force between the particles was assumed purely repulsive. 
In fact, £r ^ 1. 

The next coefficient, 



(17) 



represents the rotation of the normal vector eV during 
the collision, measured in units of ^. It is defined by 
e ^ • e ^ = cos{e^^). 

The change of the corresponding rotation velocity is 
described by a further coefficient. 



^(f) 

m 



(18) 



Using the conservation of angular momentum, L = 
^eff^^(^)^(^), we see that this coefficient is redundant 
and may be expressed through Sr- 



f(0) 



The propagation of time is accounted for by 

et = f 



(19) 



(20) 



which holds the scaled contact time. It is obviously 
needed since time is also a variable which changes during 
a mechanical contact. Its meaning becomes clear if one 
looks to the center of mass coordinate R which is not 
affected by the collision due to momentum conservation. 
To determine its post-collisional value, one needs to know 
the time when the collision terminates. 
Finally we need 



Sr 



f{f)_ 



(21) 



which is (up to the sign) the ordinary coefficient of 
normal restitution including the infiuence of centrifugal 
forces occurring for non central collisions, —£r = e^. 

Following the arguments of the previous section, the 
state of the colliding particles is completely determined 
by r(t), r(t), (f{t)^ (pit) and t. If we define 



X(0) 



(22) 
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equations (|16[)-(|21[) establish then a complete set of equa- 
tions to compute the post-collisional state, x{^)i from the 
pre-collisional one, x(0) 



We arrange the coefficients given in Eqs. (16)-(21) in 
form of the matrix of restitution 



( Sr 
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(23) 



such that the colhsion dynamics is described by the prop- 
agator 



x(t) =ex(0), 



(24) 



which has exactly the same functional form as a the tra- 
ditional propagator rule, Eq. ( 15 ) 



Similar to Eq. ( 15 ) which is the basic equation of eMD 
under the simplifying assumption, Eq. of instanta- 
neous collisions, Eq. [24] will be the basic equation of our 
generalized eMD, which does not rely in instantaneous 
collisions. 



IV. IMPROVED COLLISION RULE 

In traditional eMD simulations the particles move 
along straight lines or ballistic trajectories under the in- 
fluence of constant external fields like gravity, interrupted 
by instantaneous events (collisions) where their velocities 
are adjusted according to the collision law. That is, the 
collision law does not change the positions of the parti- 
cles. 



The new propagator, Eq. (24), requires that the cor- 
responding collision law changes both the velocities and 
also the positions of the particles. The change of the 
position in a collision will cause some problems in simu- 
lations, namely, it may happen that the designated posi- 
tions are occupied by other particles. This problem will 
be addressed in Section IVlBI 

In the present section, we detail the update of the par- 
ticles' velocities and positions, provided the new positions 
are not occupied. We describe how to apply the matrix 
of restitution e to obtain the post-collisional coordinates 
r/, V2 from the pre-collisional coordinates r2^, 

V2 for a given set of material parameters and parti- 
cle masses. For convenience, we use two (fixed) reference 
frames: The laboratory system (spanned by e^, e^, 

e^) and H as defined in Sec. |ll[ Eq. X indicates, 
that the vector X is expressed in the reference frame S. 
Vectors without a hat are expressed in S^, respectively. 



A. Position Update 

The base vectors of the laboratory frame expressed 
in S read 



(25) 



The direction of the relative coordinate after the col- 
lision reads 



= I sin(£:^^) 




(26) 



expressed in the reference frame S. The corresponding 
vector expressed in reads 



( 



6^ • 



(27) 



The distance r' between the two spheres after the colh- 
sion is given by 



r' = r°£. 



(28) 



where is its precollisional value. With this, the vec- 
tor pointing from the origin of S to particle 1 after the 
collision reads 



1712 



mi + 1712 



-r e 



(29) 



expressed in the laboratory frame S^. The corresponding 
vector pointing to particle 2 reads 



Af2' 



mi 



mi + m2 

The center of mass coordinate after the collision reads 



(30) 



R' = R° + R°etT 



(31) 



expressed in the laboratory frame. 

With this, the postcollisional particle positions ex- 
pressed in the laboratory frame read 



& + Arl 



B. Velocity Update 



(32) 



by 



The angular velocity at the instant of collision is given 



L 



meff (rO) 



(33) 
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The corresponding postcollisional value reads 

^(0) 



(34) 



The derivative of the unit vector of the postcohisional 
relative coordinate reads 



-(^' sin(£^(^) 




(35) 



in the reference frame S. The corresponding vector ex- 
pressed in the laboratory frame reads 



(36) 



V 



The normal component of the relative velocity between 
the two spheres after the collision is given by 



(37) 



where is its pre-collisional value. With this, the post- 
collisional velocity of the particles measured from the ori- 
gin of S, expressed in the laboratory frame read 



-m2 
mi 



mi + m2 



{f'e;+r'4). (38) 



With this, the post-collisional velocities expressed in the 
laboratory frame read 



vl = R'^Avl, 



In absence of external fields we have R' = R^. 



(39) 



Together with the matrix of restitution, Eq. (23), 



Equations (32) and (39) establish a complete set of equa- 



tions for the computation of the post-collisional positions 
and velocities from the pre-collisional values. 



V. COLLISION OF GRANULAR PARTICLES 

The previous section provides a general way to per- 
form event-driven simulations of soft particles, that is, 
the hard-sphere approximation, Eq. (|2| is not exploited. 
So far, however, we did not specify the particle interac- 
tion force which determines the properties of the matrix 
of restitution, Eq. (23). 



In this section we consider two widely used models for 
the interaction force F^, the linear dashpot model and 
the model of viscoelastic spheres to obtain the matrix of 
restitution, Eq. (23). Both models are characterized by 



many material and system parameters, thus, the compo- 
nents of the matrix of restitution are functions of these 
parameters. Since, for both force models, an analytical 



evaluation is not possible, by appropriate scaling we re- 
duce the problem to three independent parameters, lead- 
ing to a convenient way for computing efficient lookup 
tables for the matrix of restitution. 

Together with the collision rule, Eqs. (32) and (39), 



the results of this section allow for highly efficient event- 
driven simulation of granular gases of soft spheres. 



A. Linear-Dashpot Model 

The linear-dashpot model is widely used in the litera- 
ture for the simulation of granular systems. Its physical 
relevance may be questioned since neither the elastic [31] 
nor the dissipative part of the force [32 agree with phys- 
ical reality. It even violates a dimension analysis [18]. 
Its main characteristics is that in the hard-sphere limit it 
leads to a coefficient of restitution which is independent 
of the impact velocity (which disagrees with experiments 
as well, e.g. [33 ). Although physically questionable, 
the linear-dashpot model is widely used since its conse- 
quence, the constant coefficient of restitution, simplifies 
the analytical analysis largely. Therefore, except for very 
few examples, e.g. [34H37], virtually the entire Kinetic 
Theory of granular gases relies on this assumption. 

The linear-dashpot model defines the normal force be- 
tween colliding spheres by 



Fn = k{l -r) -jr, 



(40) 



with I = Ri -\- R2, and k and 7 being the spring constant 
and the dissipative parameter. With this force and the 
scaling (see Eq. (pT|) 



k 



(41) 



from Eq. (12) we obtain the equations of motion 



df2 



where 



X 



Cdh 



dr 
at 



(42) 



(43) 



We solve Eq. (42) with the initial conditions (see 
Eq. (|14|) 

dr 

(^(0) = 0, r{0) = l and -.(0) = -1 (44) 

dt 

for a given set of {/,C(^,Cdis} in the interval < t < f, 
where f is the time where the collision ceases given by the 
condition Eq. (fTol). The matrix of restitution, Eq. (23), is 



then obtained by using the definitions of its components, 
Eqs. ([16|-([18|, ([20]) and (l2l|. 
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The reduced set of parameters, {/, c^^, Cdis}, follows 
from both, material parameters (/c, 7, mass density p), 
particle sizes R2) and impact parameters (impact 
velocity v and eccentricity e = d/l, see Fig. |2|. 





FIG. 2. Eccentric collision of spheres. 

For practical application, we need the matrix of resti- 
tution, Eq. ([23|, for a wide range of the (physical) 
system parameters, corresponding to a certain area in 
the C(^, Cdis}-space of complicated shape. For elastic 
spheres (7 = 0) , the intervals for the physical parameters 
given in Tab. W lead to the area in the {[, c^^, Cdis}- space 
shown in Fig. Is^ showing s^^ as a function of [and c^p. 



FIG. 3. The space of physical parameters given in 
Tab. ^ translates into a space of the scaled parameters 

^Cc^/c™^^, of complex shape. The figure shows one 

of the elements of the matrix of restitution, s^, as a function 
of the scaled variables, for the special case 7 = (elastic col- 
lisions). Each point of the colored region corresponds to a 
point in the physical space given in Tab. |l] The white regions 
are inaccessible within the chosen set of physical parameters. 
Expressions for c™'"" and Z™'"' are given in Eqs. (pSl and (|43). 





unit 


min. 


max. 




k 


[10^ N/m] 


1 


1000 


spring constant 


R 


[m] 


0.001 


0.1 


particle radius 


Pm 


[kg/m^] 


250 


3250 


material density 


7 


[kg/s] 


0.01 


1.25 


dissipative parameter 


V 


[m/s] 


0.001 


25 


impact velocity 


d/l 




0.01 


0.99 


eccentricity 



TABLE I. Space of physical parameters used to obtain the 
matrix of restitution for the linear-dashpot model. For the 
definition of impact velocity and eccentricity see Fig. [2] 



with ^ —4.6. By switching to In-^, the lines of 

constant eccentricity e in the {In/, Inc^j-space are hence 
raised to straight vertical lines (see Fig. [3|. 

Further, from the definition of /, Eq. (43), the scaling 
Eq. (41) and geometry, we obtain 



luj 



(47) 



Using Eq. ( 45 ) , we express e in terms of / and and end 
up with 



We switch now from (/, c^^, Cdis) to a new set of indepen- 
dent parameters, such that the parameter space is bound 
by perpendicular straight axis. This is necessary for the 
numerically efficient access to the elements of the matrix 
of restitution (represented as a lookup table) needed for 
efficient eMD simulations. 

From the definitions Eq. ([l3), Eq. ([g]), X = //[and the 
geometry of the collision. Fig. 2| we find 



In Ceo = In r 



In 



1 



1 



(45) 



which indicates that for a given impact eccentricity, e, all 
possible {In/, Inc^j-pairs are located on a straight line of 
slope 1 [20^ with -4.6 < In - l) < 1.95 for the 

parameters given in Tab. [l[ That is, for a given [ the 
smallest accessible is given by 



In/ 



In 



In/ 



(46) 



ln[=lnf^ 



1 



In 



(48) 



For the physical parameters in Tab.|l]we obtain — 2.3 < 
In (^) < 14.83, thus, for a given impact eccentricity, the 
smallest attainable / is hence given by 



In/™ = m" 



In 



/ 



(49) 



with ^ -2.3. 

Consequently, if we would plot Fig. [3] with axis 
Inccp/c™ (instead of Inc^^) and In^r™ (instead of In/), 
the accessible data points would form a rectangular area. 
Thus, the complicated shaped colored region in Fig. [3] 
is transformed into a rectangle which allows for an effi- 
cient use of a corresponding lookup table in eMD simula- 
tions. The generalization to inelastic particles is straight- 
forward. 



7 



The inverse transformation from {In ^^jc^^^^ In^f"^^"} to 
{lnC(^,ln/} is obtained directly from the definitions of 
lnc<p/c^i- and lnyz~"^^'^: 



B. Viscoelastic Spheres 

The normal component of the interaction force be- 
tween two colliding viscoelastic spheres reads 



In / = In V m 

/mm 



f - In 



In 



(50) 



Figure [4] shows the components of the matrix of resti- 
tution, Eq. ( [23| ), for the linear dashpot interaction force, 
Eq. ( [4Q| and the range of physical parameters specified 
in Tab.^l Each row corresponds to a certain (scaled) dis- 
sipative constant Incdis, see Eq. |43) 



The first column of Fig. [4] displays (see Eq. (17)) 
which describes the rotation of the inter-particle unit vec- 
tor Cr during the contact. This rotation angel is deter- 
mined by the contact duration, r, and the rotation veloc- 
ity, (p. If we would disregard centrifugal forces, the rota- 
tion velocity would be constant. As the contact duration 
decreases with inelasticity the rotation angle and, thus, 
e^^p also decrease with inelasticity. Regarding the compo- 
nent elastic collisions hence represent the marginal 
case pD . 

The second column in Fig. [4] (see Eq. ^) shows e^. 
The coefficient en = —Sf is the well known coefficient 
of normal restitution including effects due to centrifu- 
gal forces. It describes the loss of energy of the par- 
ticles' relative velocity in normal direction, due to the 
collision. From this interpretation follows that —£f de- 
creases with increasing dissipation for all combinations 
of {In c^/c™. In [/[™}. 



The coefficient £r (see Eq. ([16])), shown in the third 
column of Fig.|4]stands for the ratio of the post- and prec- 
ollisional distance of the particles. Due to the premature 
end of collision (Eq. (10), see [16l[T7] for an in-depth dis- 
cussion) the value of Sr may differ from 1 for inelastic 
collisions (7 > 0). For impacts leading to large rotation 
velocity, the coefficient £r may significantly deviate 
from unity because of centrifugal forces. While dissipa- 
tive forces cause a premature and of the collision, they 
also reduce the rotation velocity. Consequently there is 
an optimal value for the damping coefficient, 7 (or its 
scaled value Cdis), which minimizes Sr- 

The last column of Fig. [4] shows the component St (see 
Eq. (20)) which stands for the collision duration mea- 
sured in units of the characteristic time T (see Eq. (pT])). 
The absolute value of £t decreases with damping since 
due to the premature end of collision, the contact dura- 
tion, r, decreases with increasing dissipation. 



P _ z^el I z^dis 



where 



p,,{l - rf" - -Ap,,fy/T^r , (51) 



Pel 



3(1 - 



(52) 



and y, V and i^eff denote the Young modulus, the Poisson 
ratio and the effective radius i?eff = ^1^2/(^1 + ^2), 
respectively. The elastic part of this widely used 
collision model [BHS] is given by the Hertz contact force 
[31]. The dissipative part, F^^®, was first motivated in 
j38] and then rigorously derived in ^32^ and ^ , where 
only the approach in [32] leads to an analytic expression 
for the parameter being a function of the elastic and 
viscous material parameters, see [32 for details. 

Using the normal force Eq. (51 ) and the scaling relation 
Eq. (ITT]) with 



1 



X 



(-rO) 



V5 



0' 



(53) 



where k = p/meff, the general equation of motion 
Eq. (12) reads 



dip 
dF 
d^f 
dp 



(l-r) 



3/2 



Cdis 



dr 
d! 



(54) 



where Cdis = 

Proceeding along the lines of Sec. 



V A[ we solve 



Eq. (54) with the initial conditions Eq. (44) for a given 
range of physical parameters to obtain the matrix of resti- 
tution, Eq. (23). The intervals of parameters specified in 



Table [III cover a wide range of applications. 





unit 


min. 


max. 




Y 


[10^ N/m^] 


0.01 


100 


Young's Modulus 






0.2 


0.5 


Poisson's ratio 


R 


[m] 


0.001 


0.1 


particle radius 


pm 


[kg/m^] 


250 


3250 


material density 


A 


[s] 


10-^ 


1 


dissipative parameter 


V 


[m/s] 


0.001 


25 


impact velocity 


d/l 




0.01 


0.99 


eccentricity 



TABLE II. Parameter space scanned to obtain the matrix 
of restitution for viscoelastic spheres. For the definition of 
impact velocity and eccentricity see Fig. [2] 

The set of physical parameters can be transformed in 
a set of scaled variables, {/, c^^, Cdis}- Again, the speci- 
fied ranges of physical parameters correspond to a region 




0.5 



1 



1 2 3 



FIG. 4. (color online) Components of the matrix of restitution, Eq. (23), for the linear-dashpot interaction force, Eq. (40) and 
some values of Incdis- Abscissa of all panels: In^^/c'^'''. Ordinate of all panels: Inyr"^^"^. The range of scaled parameters {c^,l) 
corresponds to the physical parameters space defined in Tab. [l] 



in the c^^, Cdis}-space of complicated shape. As in the 
case of the linear-dashpot model we look for a transfor- 
mation such that the admitted sets of parameters estab- 
lish a rectangular system. 



The parameters of the force do not enter Eq. (45), 



therefore, it holds true for viscoelastic spheres too. Since 
the marginal values for the impact eccentricity, e, remain 
(same ranges in Tabs. (l] and [n|) , again g^^^ = —4.6. 



The corresponding equation to Eq. (48) valid for the 
linear-dashpot model, reads 



In / = In 



2/5- 



In 



(55) 



for the case of viscoelastic spheres. Using the parameters 
from Tab. [ill we obtain for the first term 



0.75 ^ m™ < In 



2/5- 



< m"'^^ ^ 13.66 (56) 
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With this we define 



m 



In 



I 



(57) 



In the same way as for the hnear-dashpot force, we 



instead of 



In Coo , In / 



use [in (c^/c™), In 

as independent variables. While the domain of physical 
parameters (Tab. In]) is represented by an area of complex 



shape in the coordinates 



In Coo , In / 



(similar to 



Fig.|3|), 



in the new variables, the domain is bound by a rectan- 
gle which is much better suited for the construction of a 
lookup table for the matrix of restitution. 

In contrast to the linear dashpot model discussed in 



Sec. V A, for viscoelastic spheres the dissipative param- 



eter Cdis = f4 depends on / and c^n via T. From the 



2T 



definitions of c^^, X, T, L and geometry, we obtain 



IX c^p II c< 



dv I 

Using Eq. (|45| to replace ^/d = 
Cdis (below Eq. (54)) this yields 



^ 

vd p 



(58) 

i/e and the definition of 



Incdi, 



In 



-4 

2 / 



In/ 



1 



In 



(59) 



That is, for a given c = Ind^y), lncdis(/, c^^) forms a 
curved surface in the {/, c^^, Cdis}-space, where c ranges 
from c""^" ^ -18.71 to c""^^ ^ 9.84 for the physical pa- 
rameters given in Tab. [n| With this, we define 



^dis 



In/ 



In 



1 



(60) 



Using Incdis/c™ instead of Incdis, the physical parame- 
ters given in Tab.[ll|are mapped to a cube-shaped domain 
in the In Ccp/c;^i'^-ln yr"^^"-ln cdis/c™-space, allowing for effi- 
cient lookup tables. 

Fig. [5] displays the result for a selection of dissipa- 
tive parameters Incdis/c™. Similar to Fig. |4] each row 
of Fig. [5] shows the four components of the collision map- 
ping Eq. ( [23| ) for a fixed dissipative parameter In cdis/c^^i^^. 
Again, dissipation increases from the top to the bottom 
row. The discussion of Fig. [5]is absolutely equivalent to 
the linear dashpot case, FigT^ 

The corresponding transformation back to /, and 
Cdis may be obtained directly from the definitions and 
reads 



In [ = In ^J— 

/mm 



- m 



In Ccr 



In Cdis = In 



ln/ + ^/" 



■-dis 



In 



In [ In 



2 ln3^SF+^" 



1 + e 



2 ln-^+^" 



(61) 



VI. EVENT-DRIVEN MOLECULAR 
DYNAMICS ALGORITHM 

A. Traditional event-driven Molecular Dynamics 

The traditional eMD scheme of hard particles is rather 
simple although an efficient implementation allowing for 
the simulation of many millions of particles may be tech- 
nically rather complex, see e.g. [15^. Its basic concept is 
to 

1. ffnd the next colliding pair (i, j) of particles in the 
system and their collision time t* 

2. propagate all particles k to this time, 

rk:=fk^Vk{t"-t) (62) 

where t is the present time and is the present 
velocity of particle k 

3. compute the post-collisional velocities of particles i 
and j due to the collision rule 



(63) 



where Sn is the coefficient of normal restitution. For 
simplicity of the notation we consider here particles 
of identical mass, the generalization is straightfor- 
ward. 

4. continue with step 1 

While all eMD schemes work in principle as described, 
there are many ways to increase the efficiency, e.g. 
[El [14[ [151 Ho] , which shah not be discussed here. More- 
over, the scheme described above does not take into ac- 
count external fields like gravity, interaction with (mov- 
ing) boundaries etc. 

Figure [6] (top) shows the trajectories of two colliding 
particles as obtained by eMD in comparison to force- 
based MD, that is, the numerical solution of Newton's 
equation of motion. For the interaction force we as- 



sume a linear-dashpot model, Eq. (40), with the pa- 
rameters k = 2 kN/m, R = 0.1 m, p = 1140 kg/m^, 
'U = 5m/5, 6 = 0.3 (see Fig. Initial velocities are 
= ('^/4, v/2^0) and V2 = ('^/4, —v/2^0). Since here 
we assume elastic interaction (7 = 0), the corresponding 
coefficient of normal restitution is = 1. 

The figure reveals two fundamental problems which 
are both attributed to the assumption of instantaneous 
collisions: First, the finite duration, r, of collisions in 
the physical system leads to a finite rotation of the inter- 
particle unit vector. Consequently, the directions of the 
final velocities differ for MD (based on forces) and eMD 
(based on the coefficient of restitution). As indicated in 
the figure, the deviation may be large. Only for the case 
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FIG. 5. Components of the collision mapping Eq. (23) for the viscoelastic interaction Eq. (51 ). Abscissa of all panels: \n^<p/c 
Ordinate of all panels: In^f" 



^. Each row displays s^p, Sr, Sr and et for the dissipative parameter In^dis/c 
white label in the corresponding image for e^. Parameters as indicated in Tab. [Tl| (color online). 



indicated by the 
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FIG. 6. (color online) Traces of two colliding spheres (param- 
eters are given in the text). Black lines show the numerical 
integration of Newton's equation (MD), red lines (top) show 
the trajectories as obtained from eMD with the assumption 
of instantaneous collisions. The green lines (bottom) show 
the trajec tories as obtained by the new eMD algorithm (see 
Sec. VI B). Symbols and numbers (of the respective color) in- 
dicate the particle positions at equidistant points in time. The 
number stands for the moment when the particles touch and 
7 corresponds to the end of the collision (stepsize dt = t/7). 
The dashed circles show the spheres at the moment of impact. 



of a central collision, the directions of the final velocities 
agree for MD and eMD. 

Second, the position of the particles as a function of 
time may be different for MD and eMD. This applies 
to both, central and off-central collisions. In Fig. [6] we 
indicate the dynamical properties by plotting dots (of 
the respective color) on top of the trajectories (lines) at 



equidistant intervals of time. 

We wish to mention that the chosen parameters for the 
plot in Fig. |6] correspond to very soft particles in order to 
visualize the differences between MD and eMD. A care- 
ful analysis [20, 21 shows that the differences may be 
large also for more realistic material and system prop- 
erties. The fundamental problems detailed above are 
always present when collisions of physical particles are 
modeled by eMD assuming instantaneous collisions. 



B. Improved event-driven Molecular Dynamics 

i. Classification of events 



The propagation rule, Eq. (63), is used in traditional 
eMD and relies on the coefficient of restitution, e^, and 
instantaneous collisions. It shall now be replaced by 



the propagation rule Eq. (24) using the matrix of resti- 
tution, £, which takes the finite duration of collisions 
into account. We propose an improved eMD algorithm 
which does not suffer from the problems described above, 
caused by the assumption of instantaneous collisions. 

In the improved eMD scheme (eMD*), each collision is 
represented by 3 instantaneous events. These events may 
be of type E'^ or E^: Assume two particles (i^j) collide 
at time t*. This collision is represented by 

a) an event of type E'^ at time t* where the positions 
of the particles are set due to the rotation of the 
inter-particle unit vector. The velocities are set to 
the center-of-mass velocity of the colliders, 

b) an event of type E^ at time + r where particle i 
adopts its post-collisional velocity, and 

c) an event of type E^ at time t* + r where particle j 
adopts its post-collisional velocity. 



2. Events of type 

An event of type E"^ occurs at the moment t* when 
a pair of particles (i,j) gets in contact, similar to the 
events in traditional eMD. The following sub-tasks are 
performed: 

1. Compute the scaled parameters (In/, Inc^^, Incdis) 
from the physical material parameters, the particle 
radii, the impact geometry and the velocities of the 
particles. 

2. Compute the components {e^p^St^Sr^ Sr} of the ma- 
trix of restitution, Eq. (23). This may be done in 



a convenient and efficient way using lookup tables 
based on the transformations described in Sec. [Vl 

3. Apply the collision rule detailed in Sec. |IV| with 
= to rotate the particles around their center 
of mass by the angle ip = e^^^ and to obtain the 
postcollisional velocities, vl and v-. 
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4. Set the velocities of both particles to the center of 
mass velocity Vi = Vj = R^. 

5. Store the computed post-collisional velocity i// in a 
local variable of particle i and, respectively, Vj in a 
local variable of particle j. 

6. Mark both particles as collision not yet accom- 
plished by setting a local flag. 

7. Schedule two more events of type in the global 
event list, both occurring at time +r; (r = £tT): 

(a) the velocity of particle i will be updated. 

(b) the velocity of particle j will be updated. 



3. Events of type 

If a particle suffers an event of type E"^ at time t*, it 
suffers an event of type E^ at time + r, where r is the 
duration of the collision which was computed when the 
event of type E"^ was handled. In difference to events 
of type E'^ describing two-particle interactions, events of 
type E^ concern only one particle. The following sub- 
tasks are performed when an event of type E^ occurs: 

1. Check whether the flag collision accomplished is set 
in the concerned particle i. If this is the case, do 
nothing. Otherwise continue with item 2. 

2. Set the velocity of the concerned particle i to the 
value which was previously computed and stored 
in a local variable of particle i, see item [5] in Sec. 
IVIB21 

3. Set the flag collision accomplished in particle i. 



2. if t* is smaller than the first (next in time) entry 
in the collision list, propagate all particles k to this 
time. 



rk :=rk^Vk{t'' - t) , 



(64) 



handle the collision as an event of type E"^ and 
proceed with step 1. 

3. propagate all particles k to the time of the next 
scheduled event of type E^ 



Tk := Tk ^Vk{t^ - t) , 



(65) 



handle this event and remove the entry from the 
list. If there is more than one event scheduled for 
the same time, chose any of them. Proceed with 
step 1. 

An exemplary application of the algorithm is shown in 
Fig. [6] (bottom). The black lines again denote the trajec- 
tories due to Newton's equation (same as upper panel). 
The green lines display the trajectories as obtained by 
the eMD* algorithm. At time when the particles get 
in contact, an event of type E"^ is performed. This event 
rotates the inter-center unit vector around their center of 
mass, and, thus, relocates the particles instantaneously 
to new positions (time is shown twice). From there on, 
the particles move at the velocity of the center of mass. 
At time 7 two events of type E^ occur where both parti- 
cles adopt their final post-collisional velocities. From this 
time on the trajectories due to eMD* and MD (Newton's 
equations) agree perfectly. In contrast, as indicated by 
the upper panel of Fig. [6] the results of eMD (red lines) 
and MD differ significantly. 



5. Exceptions 



4' Schedule of events 

The eMD* algorithm is similar to the eMD algo- 
rithm in the sense that the computation proceeds from 
one event to the next. The particle velocities are only 
changed due to these instantaneous events (except for the 
trivial acceleration resulting from homogeneous and con- 
stant external fields which does not influence the collision 
sequence). In eMD*, the events of type E"^ correspond 
to the events in eMD. 

Again, we discuss only the principle of the algorithm, 
not the technicalities of its implementation. We assume 
there is a global list which contains the sequence of sched- 
uled events of type E^. Initially, the list is empty. 

The eMD* algorithm then works as follows: 

1. find the next colliding pair (i, j) of particles in the 
system and their collision time t* (begin of the col- 
lision) 



For the description of the algorithm, we silently as- 
sumed that the operations due to events of types E'^ and 
E^ are permitted. This is, however, not always the case 
but we have to deal with two possible exceptions: 

a) The ro tation step (event of type E'^, item [s] of 



VIB2) may not be executed as it would lead to 



overlap with other particles. 

b) In the time interval between an event of type E'^ 
and the associated events of type E^ the colliding 
particles (i, j) move at the center of mass velocity, 
which is an unphysical but very short-lived tran- 
sient state. In this time interval, one of the particles 
(i, j) (or both) may collide with another particle k. 

Case a): We are of the opinion that event-driven MD 
is restricted to the domain of dilute systems. To apply 
eMD, we have to assure that in the corresponding phys- 
ical (force-based) system, the frequency of three-particle 
interactions is negligible as compared with the frequency 
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of two-particle interactions (see \W for a detailed discus- 
sion of this problem). For realistic and relevant material 
and system parameters, the rotation angle of the inter- 
center unit vector, cp^ is rather small. Consequently, only 
minimal extra space is needed to perform this rotation. 
The probability that this small rotation would lead to 
overlap with other particles is, hence, small as well. 

In case that such an exception occurs, we fall back to 
the traditional eMD scheme for this particular collision: 
Only the velocities of the particles are changed accord- 
ing to the collision rule, Eq. (63), using the coefficient of 
restitution En = —Sr but not the full matrix of restitu- 
tion. 

Case b): Assume the collision (i, j) at time t = re- 
quires for both particles i and j an event of type at 
time r. Assume further that another particle k collides 
with i at time < r. In this case we exceptionally per- 
form the event of type E^ of particle i at time t^, just 
before the event of type E"^ of the pair (z, /c). 

That is, only the instant in time when the events of 
type E^ are executed is modified. The post-collisional 
velocities are not affected by the exception handling and, 
hence, neither conservation of momentum nor conserva- 
tion of energy are violated by this type of exception. But, 
as a consequence of the exception, the events of type E^, 
both scheduled at time r originally, are no longer executed 
simultaneously due to the interference of a third particle. 
This, in turn, violates conservation of angular momen- 
tum by a tiny amount. However, it may be shown, that 
any application of periodic boundary conditions leads to 
much stronger violations of angular momentum. 

The tiny violation of the conservation of angular mo- 
mentum may actually be avoided: The events E^ of both 
particles i and j, originally scheduled at time r are ex- 
ecuted at the earlier time t^ when either particle i or 
j interferes with a third particle k. However, this re- 
quires post-collisional communication between the par- 
ticles i and j being an algorithmic complication which 
may cause significant loss of computational performance. 
Depending on case specific demands, one has to decide 
between absolute accuracy and maximal algorithmic ef- 
ficiency. But anyway, for dilute systems, being the scope 
of the eMD* algorithm, eMD* including post-collisional 
communication is still by orders of magnitude more effi- 
cient than force-based MD. 

Furthermore, the exceptions of type b) implicate the 
question of how to deal with collisions being interfered 
by more than a third particle. For a gas considered here, 
the mean free time is much larger than the time lag be- 
tween the events of type E"^ and E^, corresponding to 
the duration of a collision. Therefore, if the frequency of 
an exception is small (~ 0.1%, see Sec. VIB6), the prob- 



6. Confidence Regions of the eMD* algorithm 

Aim of the eMD* algorithm is to simulate soft spheres 
while maintaining the advantages of event-driven model- 
ing which, in its traditional form, relies on hard sphere 
interaction. Of course, this goal may only be achieved 
if the unavoidable exceptions detailed in | VI B 5 



and hence negligible events. In this section we, therefore, 
assess the range of validity of the eMD* algorithm by pro- 
viding statistics on the exception frequency. To this, we 
simulate a granular gas of N = 10 000 elastic particles 
(interaction-force Eq. (51), A = s, pm = 1140 kg/m^, 
= 0.1 m, = 0.4). As simulation setup we choose a 
periodic box of Volume V^sim, in which the particles are 
initially located on a crystal lattice, from which they are 
then released to move freely with random velocities dis- 
tributed in a way that the resulting thermal velocity is 
about 2 m/s. 

Obviously, the frequency of both exceptions (type a, 
rotation step impos sible, a nd b, three particle contact, 
see enumeration in 



packing fraction 



VI B 5 ) is mainly governed by the 



3 V^im 



(66) 



ability of a four-particle interaction is even much smaller 
(~ 10~^). Hence, these cases may safely be neglected. 



and the Young's modulus of the particle material, which, 
in turn, infiuences on the contact duration and the rota- 
tion angle (p = e^p^^ respectively. During the simulation 
we record the number of exceptions of type a) and b) for 
various packing fractions and Young's modulus ranging 
from very soft materials like e.g. rubber to hard mate- 
rials like e.g. glass. The result is shown Fig. [?[ First 
we see that in the limit of very hard spheres or very di- 
lute systems the probability of both exception types van- 
ish. Second, for system parameters typically used in the 
literature on granular gases and for common materials, 
exceptions of both types are rare events (about 0.1%). 

Clearly, the percentage of collisions where the eMD* 
algorithm fails (and we fall back to the traditional colli- 
sion rule) is small. The eMD* algorithm, hence, indeed 
improves trajectory accuracy for typical systems. 



VII. SUMMARY 

Basic concept of event-driven Molecular Dynamics is 
the assumption of perfectly hard spheres leading to in- 
stantaneous collisions, such that the particle positions 
do not change during the collision. This assumption al- 
lows to describe the dynamics of a granular system as 
a series of independent binary collisions. Each of these 
collisions is modeled by a simple multiplication of the 
pre-collisional relative velocity in normal direction with 
the coefficient of restitution to obtain the post-collisional 
value and finally the post-collisional vectorial velocities. 
The only parameter characterizing the collision is the co- 
efficient of restitution, containing all the physics of the 
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0.01 
rubber 



Y (GPa) 



FIG. 7. (color online) Exception frequency fe in i/iooo for a 
free granular gas of elastic particles as a function of a) the 
packing fraction 77 at y = 1 GPa (upper panel) and b) the 
particles Young's modulus Y at 77 = 0.001 (lower panel). See 
text for details on the setup and the parameters. 



particle interaction. For central collisions, it can be de- 
rived by integrating Newton's equation of motion for an 
isolated pair of colliding particles (which may lead to a 
velocity dependent coefficient of restitution), e.g. [16]- 
[T9] . The coefficient of restitution is then found from its 
definition, Eq. ([3|. Hence, for central collisions, an event- 
driven description yields the correct post-collisional ve- 
locities if compared to the integration of Newton's equa- 
tion of motion. However, even for central collisions, the 
temporal properties are not correctly reproduced since 
the finite duration of collisions is neglected within event- 
driven modeling. 

Clearly, the assumption of instantaneous collisions is 
an approximation: Physically, the trajectories are deter- 
mined by Newton's equation of motion with appropri- 
ate forces and material parameters. Any instantaneous 
change of the velocities would correspond to diverging re- 
pulsive forces between the particles. Furthermore, the re- 
quest for a loss of energy of colliding particles (expressed 
by the coefficient of restitution) is not consistent with the 
assumption of instantaneous collisions since otherwise a 
finite amount of energy must be dissipated in vanishing 
time. That is, the hard-sphere model may be inappro- 



priate for the description of dissipative systems. 

The finite duration of physical collisions leads always 
to a finite rotation of the inter-particle unit vector e^. 
Only for central collisions eV remains unchanged. For the 
case of adhesive nano-particles [41 it was recently shown, 
that at very large impact rate the rotation of eV may be 
large. This rotation, in turn, causes a large deviation be- 
tween the trajectories as obtained in eMD (applying the 
collision rule, (|63|) and MD (integrating Newton's equa- 
tion). This result was generalized to oblique collisions of 
particles interacting via any force law [20 , 21 . 

Consequently, due to the hard- sphere assumption, 
eMD agrees with MD neither regarding the spatial nor 
the temporal properties of the trajectories. The devia- 
tions may be large |20l [21] . 

In the present paper we propose an alternative event- 
driven algorithm, eMD*. The essence of the eMD* al- 
gorithm is an extended collision rule. In contrast to the 
one of classical eMD, it changes not only the particle ve- 
locities but also their positions. Pre- and post-collisional 
states of the system differ in more than just the normal 
component of the relative velocity. We arrange all chang- 
ing quantities in a vector which completely describes the 
system state. Compared to classical eMD, where pre- and 
post-collisional normal component of the relative veloc- 
ity are related by the coefficient of restitution, pre- and 
post-collisional state vectors are, consequently, related by 
a matrix within eMD*. We termed this matrix matrix of 
restitution. Together with the concept of the state vec- 
tors, it allows to maintain the mathematical form of the 
hard sphere collision rule applied within classical eMD. 
Similar to the coefficient of restitution in eMD, all phys- 
ical properties of the collision are mapped to the matrix 
of restitution. The eMD* algorithm does not assume in- 
stantaneous collisions. If applicable, the post-collisional 
particle positions and velocities obtained by eMD* agree 
with those obtained by integrating Newton's equations. 
Algorithmically, in eMD* each collision is represented by 
3 events of two different types which together map the 
pre-collisional state to the post-collisional one. 

Centerpiece of the method is the setup of the matrix 
of restitution as a functional of the particle interaction 
force law. We apply the eMD* algorithm to two ex- 
amples which are important for practical applications, 
the linear-dashpot force and the viscoelastic Hertz force. 
Both force laws are characterized by two material prop- 
erties. The geometry of the particles and the vectorial 
pre-collisional velocities are further parameters describ- 
ing the impact. For both examples we demonstrate that 
the collision can be fully described by a set of three pa- 
rameters which allows to represent the elements of the 
matrix of restitution in the form of universal lookup ta- 
bles. Using these tables the eMD* algorithm turns into 
a very efficient simulation method. 

We applied the eMD* algorithm to the oblique colli- 
sion of two spheres and obtain identical post-collisional 
velocities as compared with Newton's equations. The 
trajectories are identical as well, except for a short-lived 
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transient state whose duration is of the order of the dura- 
tion of the cohision. This means, that for dilute systems, 
where the exceptions detailed in |VI B 5| are rare, negli- 



gible events, eMD* simulations are equivalent with MD 
simulations. In fact, as shown in VIB 6[ the frequency 
of (algorithmic) failure may be reduced to any desired 
number by reducing the system density, while the physi- 
cal effects of finite interaction forces are preserved. That 
is, both methods simulate granular systems composed of 
soft spheres and yield the same trajectories as functions 
of time. At the same time, as an event-driven algorithm, 
eMD* is much more efficient than force-based MD. So 
far, we only considered frictionless interactions. This, 
however, is not a principal restriction and extending our 
findings to rough, frictional spheres is subject of future 
investigation. 

Besides standard eMD, also the Kinetic Theory of 
granular gases is based on the hard sphere model since 
the Boltzmann equation is applicable only for hard 
spheres. This raises the question how the deviations be- 
tween the trajectories obtained by means of the coeffi- 



cient of restitution and from Newton's equation, affect 
the results of Kinetic Theory like, e.g., transport coeffi- 
cients, which is subject of current research. For granular 
gases it is known that the vectorial particle velocities 
are correlated due to the dissipative nature of the inter- 
actions which necessarily implies a violation of molecu- 
lar chaos [42H45] . It may hence be expected that the 
improved trajectory accuracy achieved by the eMD* al- 
gorithm is not screened by Molecular Chaos and leaves 
its fingerprints also in measurable macroscopic quantities 
like e.g. the coefficient of (self-)diffusion. 
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